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r- 1 In this paper, based on the homotopy continuation method and the interval 

Newton method, an efficient algorithm is introduced to isolate the real roots 
of semi-algebraic system. Tests on some random examples and a variety 
of problems including transcendental functions arising in many applications 

j— i show that the new algorithm reduces the cost substantially compared with 

^ the traditional symbolic approaches. 
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1. Introduction 

m > 

The problem of counting and isolating real solutions of nonlinear system 
is an important topic in computing geometry and many other applications in 
various fields, e.g., the real intersection points for piecewise algebraic curve 
PQ El El IH E], the stability of a large class of biological networks [5J [7J, 
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discovering non-linear ranking functions of loop programs in [SJ[H], automated 
proving inequality type theorem pH] and so on. 

Many algorithms for real root isolation of one polynomial in one variable 
have been developed in [TH [121 1131 Ell US] • For multivariable case, Xia et 
al.[TB] proposed an algorithm based on Wu's method for isolating the real 
roots of semi-algebraic system with integer coefficients, and made it more 
available with interval algorithm in their later work [17] . There are also 
other algorithms based on different techniques, see [HI [191 EDI EE] for more 
details. 

Actually, most of the algorithms mentioned above can compute the exact 
results because they depend on symbolic computations, but they are re- 
stricted to small size systems because of the high complexity of the symbolic 
computation. In order to avoid this problem, Shen et al. [22] presented a nu- 
merical algorithm improving the efficience based on homotopy continuation 
method combined with interval Newton iteration technique. 

In this paper, we extend the numerical method to a class of semi-algebraic 
systems and transcendental functions. We always assume that the system is 
square and it only has isolated roots in C. In order to avoid the singularity 
of the Jacobian matrix, we also suppose that the multiplicity of these points 
are one. 

The rest of the paper is organized as follows. Definitions and preliminaries 
about homotopy continuation method and interval arithmetic are given in 
section 2. Section 3 first gives a brief review of the numerical algorithm for 
isolating real roots of a class of polynomial systems, and then presents the 
numerical method for isolating the real roots of semi-algebraic systems. The 
experimental results together with comparison to symbolic methos is given 
in this section. Applying the new algorithm to some problems arising from 
piecewise algebraic curve, chemical engineering, robot kinematic problem, 
circuit design are shown in section 4. Finally, section 5 draws a conclusion 
of this paper. 

2. Interval Arithmetic And Homotopy Continuation Method 

In this section, some basic theories and tools about interval arithmetic 
and homotopy continuation method are presented. 
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2.1. Interval Newton's method 

Interval operations were first introduced by Moore [23]. The key idea of 
interval arithmetic is substituting an interval for a floating number in numer- 
ical computation, which is used to tackle the instability and error analysis. 
For a more detailed and complete discussion the reader is referred to [21] . 

A real interval X is a nonempty set of real numbers 



the set of all intervals in R is denoted by 1(E). 

An interval vector X is a vector whose elements are intervals, and an 
interval matrix can be similarly defined. In the rest of this paper, we let X 
denotes the set of interval vector. 

For an interval X, the midpoint of X is m(X) = (x + x)/2, the width 
of X is w(x) = x — x, and the radius of X is r(X) = (x — x)/2. 

Given interval X — [x, x], Y — [y,y], the four element operations are 
defined as 

X + Y = [x + y, x + y] 
X - Y = [x - y, x - y] 

X ■ Y = [min{xy, xy, xy, xy}, max{xy, xy, xy, xy}} 
X/Y=[x,x]-[l/y,l/$,Ot[y,y]. 

The intersection of two intervals X and Y is empty if either x < y or 
y < x. In this case , we write X R Y = 0. 

For interval matrices and interval vectors, the concepts such as midpoint, 
width, radius, etc, and the arithmetic operations are defined in components. 

Let R[je] be the ring of polynomials in the variables (xi, X2, • • • , x n ) with 
coefficients in M, and /= [fx, ■ ■ ■ , f n ] be a polynomial system, where f\ G 
R[x]. 

Suppose / : W 1 — > R is a function from a real vector to a real number, 
F : J(IR n ) is called an interval extension of / if 



X = [x, x] = {x G R : x < x < x}, 



F([x!,Xi], ■ ■ ■ 




for all X{ G Xi, i — 1, 2, . . . , n. 

We say that F = F(Xi, ■ ■ ■ ,X n ) is inclusion monotonicity if 



YiCXi,i 



1, • ■ ■ , n =^ F(Y 1 , ■ ■ ■ ,Y n )CF(X 1 ,--- ,X n ). 
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And it is easy to prove that all the polynomial operations satisfy the inclusive 
monotonicity. 

In the following, we let /' be the Jacobian matrix of / , and F, F' be the 
interval extension of / and /' with inclusive monotonicity, respectively. 
Moore first define the following interval Newton's operators: 

N(X)=m(X)-V(X)f(m(X)), (1) 

where V(X) is an interval matrix containing .F"(X) _1 . 

In order to avoid the computation of the inversion of interval matrix in 
formula (1), Krawcayk proposed the following operator: 

K(y, X) = y - Y/(y) + (I — YF '(X))(X - y), (2) 

where y is chosen from the region X, I denotes the unit matrix, and Y is an 
arbitrary nonsingular matrix. 

It has been proved that Krawcayk operator has the following properties. 

Proposition 2.1. Suppose K(y,X) is computed by formula (2), then 

1) : x* E X and f (af ) = => x* G K(y, X). 

2) : K(y, X) C int(X) =>■ f has only one root in X, where int(X) denotes 
the topological interior of the box of X. 

3) : K(y, X) C int(X) =>■ / has a root in X. 

4) : K(y, I)nl=04 no solution in X. 

In particular, if y and Y are chosen to be y = m(X) and Y = [m(F' (X.))] -1 
respectively, then the Moore form of the Krawcayk operators is : 

K(X) = m(X) - [m{F / (X))]- 1 /(m(X)) + A, (3) 

where A = (I — [m(F '(X))]-^ '(X))(X - m(X)). 

2.2. Homotopy continuation method 

Homotopy continuation method is an efficient numerical method for find- 
ing all isolated solutions of polynomial system. The method traces a path 
from the solution of an easy problem to the solution of the given problem by 
use of a homotopy continuous transformation . See reference [261 EH EHJ [291 
130] for more details. There also exists some software packages [311 [321 [33] for 
homotopy continuation methods. In our implementation, we use Hom4ps-2.0 
which could return all the approximate complex zeros of a given polynomial 
system efficiently, along with residues and condition numbers. 
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3. Algorithm for isolating real zeros of semi- algebraic system 

In this section, we will present an algorithm for isolating the real roots 
of a zero-dimensional semi-algebraic system. Our idea is to compute the 
isolated real root intervals of the zero- dimensional polynomial system first. 
In this step, we will use the hybrid algorithm introduced by Shen et, al [22J. 
Then decide whether these intervals satisfy inequations by substitution. 

First, we introduce the definition of semi- algebraic system. A semi- 
algebra system denoted by SAS can be written in the following form: 

A = o,.-. ,/ n = o, 

Pi > 0, ••• ,p q > 0, () 

ri! > (),-•• ,n s > 0, 1 ' 

where n > 1 and q, s, t > 0. We call it zero-dimensional semi-algebraic system 
if {/i, • • • , f n } has only finite zeros in C. Following the notations in software 
package Discover [16], we let /, n, p, h denote the polynomial equations, 
non-negative polynomial inequalities, positive polynomial inequalities and 
polynomial inequations respectively. 

3.1. Real root isolation for zero- dimensional polynomial system 

In this subsection, we introduce a hybrid method for the real root isolation 
of a zero- dimensional polynomial system, see [22] for more details. 

Suppose we have obtained all the isolated approximate roots of / through 
homotopy continuation method. For each approximate root, the following 
proposition gives a method to construct initial interval which contains its 
corresponding accurate root. 

Proposition 3.1. [22] Let f = ■ ■ ■ , f n ] be a polynomial system in M.[x], 
and x e R n be an approximate zero of f. If the following conditions hold: 
1: f'(x) exists, and there are real numbers B and r] such that 

2: There exists a ball neighbourhood 0(x,u) such that f'(x) satisfies the 
Lipschitz condition on it: 

11/ '(as) -f'(y)\\ < K\\x-y\\yx,yeO(x,u), 
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3: 

1 1 - 2^1^2h 
h = BKu < - , u > u, 



then f has only one root in 0(x,u). 

Let x be an accurate root, and x be its approximation. Denote the 
Hessian matrix of fj by Hj, and h % - its column vector. Let 

X = max 1<t<n y\hUx)\ max , andr= 11 ^^Hl^^n > ( 5 ) 

where | • | max denotes the maximum module component of a vector. 

A natural fact about r is that if r > Im(x), then s is a imaginary number, 
where Im(x) means the imaginary part of x. So we can delete some complex 
roots without using the interval arithmetic. 

We use algorithms from [22] to compute u and A. 

According to the above descriptions, the following algorithm real -roots 
can be used for isolating real root intervals of /. 
Algorithm 1 : real -root -isolate [22] 

Input: Polynomial system: /, and a threshold r. 

Output: Isolated intervals of / real -roots or {}. 

1. Let real _roots — {}. 

2. Computing all the isolated roots of /, denote it by roots, and let ureal 
be the number of roots. 

3. If ureal = 0, then return real-roots, and stop this program. 

4. For 2 = 1: ureal; 

• let z be the zth element of roots, and determine whether z is a 
complex root using formula (5) and Im(z). 

• If z is a complex root, stop; else, construct the initial interval of 
the real part of z according to proposition 3.1, and denote it by X . 
Compute K(X.o) according to formula (3). If K(X ) C X , then 
let real-roots = real-roots U {K(~X ) fl X }, and stop this step. 
If K(X ) H X = 0, stop. Otherwise, bisect X = K(K) n X , 
and process each half separately. 

5. End For 
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6. If there exists some intervals in real-root such that they are not dis- 
joint, then compute the intersection of these intervals, split each inter- 
val Xj = Xji UXj2, where Xji denotes the intersection portion, process 
all the subintervals, and remove the subinterval which don't contain a 
real root of /. 

7. If the width of an interval is bigger than r, then bisect the interval and 
process each half separately. 

8. Return real-roots. 

3.2. Algorithm for semi- algebraic system 

Now suppose we have obtained the isolated real root intervals of polyno- 
mial system / defined in (4). Let X = ([ai, &i], • • • , [a n , b n ]) be an isolated 
real root interval of /, and x be the accurate zero of / such that x G X. 
Assume / is a positive polynomial inequality in (4) and /(X) = [a, b], then 
there are only three cases that can happen about the sign of [a, b], including 
the following: 

1: a > f(x) > 0. 

2: b < =>• f(x) < 0. 

3: G [a, b]. 

Hence, if [a, b] satisfies the first case or the second, we can easily decide 
whether X satisfy / or not. For case 3, it is difficult to decide the sign of 
f{x). In the following, we present a method to solve this problem. 

Let g = y 2 f + 1 be a new polynomial in R[cc,y], and [f,g] be a poly- 
nomial system {/i, • • • ,f n ,g} in K[cc,y], where R[a;,y] denotes the ring of 
polynomials in the variables (x±,X2, • • • , x n , y) with coefficients in R. Assume 
Z = {Zi, • • • , Z m } is the set of the isolated real root intervals of [/, g\. Define 
the projection 7r : J(IR n+1 ) — > 7(R n ) which remove the last coordinate of an 
interval vector in 

Theorem 3.1. If the intersection of X and n(Zi) is an empty set for all 
% G {1, 2, • • • , m}, then f(x) > 0. 

Proof. Denote X by ([ai, h], • • • , [a n , b n \) and Z { = ([an, bn], ■ ■ ■ , [a in , b in ], 

[yn,yi2]),i = 1,2, • • • ,m. 

Suppose f(x) < 0, then it is easy to see that the point (x, yj—l/f(x)) 
is the accurate zero of [/,<?], so there exists an isolated real root interval Zj 
which contains the point (x, y/—l/f(x)). Hence the intersection of X and 
7r(Zj) is not an empty set, which contradicts to our suppose, so f(x) > 0. 
This proves the first part of the theorem. 
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On the other hand, assume f(x) > 0, then g = y 2 f(x) + 1 > 1 for 
any real value of y, so n(Zi) does't contain x for all % e {1, • • • ,m}. This 
establishes that the intersection of X and -n(Z.j)(i — 1, • • • , m) is an empty 
set. This completes the proof. □ 

Similar to theorem 3.1, if we construct g = y 2 f — 1, then we can decide 
whether f(x) < or not. 

According to the above theorem, the following Deter -Sign can be used 
to determine the sign of f(x) if G [a, b]. 
Algorithm 2: Deter _Sign 

Input: A polynomial system /; An isolated real root interval X of /; A 
polynomial /, and an index /, the value of / is 1 or -1. 

Output: If f(x)I > 0, then return 1; else return -1; 

1. If 1=1 

• Let g = y 2 f + 1; 

2. Else 

• Let g = y 2 f - 1; 

3. End If 

4. real^roots = real_root^isolate{\f , g]); 

5. sum = 0; 

6. For i = 1 : lengthireal -roots) 

• If X f] n (real -roots {i}) = 0, then sum = sum + 1; 

• Else, return -1; 

• End If 

7. End For 

8. If sum = ureal then return 1; 

9. End If 

Using algorithm 2, it is easy to see that if 

Deter _sign(f , /, X, 1) = 1 A Deter-sign(f , f,X,—l) = 1, 
then f(x) = 0. 

In the following, a simple example is given to explain how the algorithm 
determines the sign of f(x). 
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Example 3.1. Let f— [x 2 + y — 2, x + 2y — 3], and h = 3x + y — 4. /Hs easy 
to see that f and h has a common point (1, 1). 

Using algorithm 1, we obtain an isolated real root interval of / which 
contains the point (1, 1). 

X 1 = [0.999999999999900,1.000000000000100], 
[0.999999999999900,1.000000000000100]. 

Substituting X\ into h, we have 

= 10 -11 x [-0.400124378074906,0.399680288865056]. 

Now, let hi = a 2 (3x + y — 4) + 1, then the isolated real root intervals of 
[/,/ii] are 

Zi = [ - 0.500000000000100, -0.499999999999900], 
[1.749999999999900,1.750000000000100], 
[ - 0.516397779494422, -0.516397779494222]. 

Z 2 = [ - 0.500000000000100, -0.499999999999900], 
[ 1 . 749999999999900 , 1 . 750000000000 100] , 
[0.516397779494422,0.516397779494222]. 

So Xi n 7r(Zi) = and X 2 H tt(Z 2 ) = 0, then we have h(l, 1) > 0. 
Similarly, we obtain h(l,l) < through isolating the real roots of [ / , a 2 (3x+ 
y-A) -1]. Hence /i(l,l) =0. 

Now, we give the following algorithm for computing the isolated intervals 
which satisfy the inequations h in (4). 
Algorithm 3: Dele-Inequ 

Input: Polynomial system / and h; isolated real root intervals real-roots 
of/. 

Output: real-roots = {X £ real-roots\h(x) ^ 0,f(x) = 0, x £ X}. 

1. Set Index = {}; 

2. For s = 1 : length(real-roots) 

• Substituting rea/_ roots {s} into h, and suppose 
/i(rea/_roots{s}) = ([ai, &i], • • • , [a t , 6 t ]); 

• Let set = {i £ {1, • • • , t}|0 £ [a,, 6i]}; 

• If length(set) ^ t, then break; 
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• End if 

• For t — 1 : length(set) 

— Sigrii = Deter ^sign(f ,real~roots{s}, h t ,l); 

— If sigrii = 1, then 

Sign 2 = Deter _sign(f ,real-roots{s}, h t , —1); 

* If sigri2 = 1, then Index = Index {J{s}, break; 

* End If 

— End If 

• End For 

3. End For; 

4. real^roots = real-roots \ real -roots {Index). 

In the following, we will present a method to deal with the non-negative 
systems n in semi-algebraic system defined in (4). 
Algorithm 4 : Dele-Nonnega 

Input: Two Polynomial systems: / and n; isolated real root intervals of 
/: real-roots . 

Output: real_roots = {X £ real_roots\n(x) > 0,f(x) = 0, x G X}. 

1. Set Index = {}; 

2. For t — 1 : length(real_roots) 

• Substituting real-roots{t} into n, and suppose 
n(real-roots{t}) = ([ai, h], ■ ■ ■ , [a s , b s }); 

• If there exists some i G {1, ••• ,s} such that 6j < 0, then let 
Index = Index U {t}; 

• Else if there exists some % G {1, • • • , s} such that G [a,, Let 
set = {ie{l,-- - , s}|0 G [a^ &»]}; 

— For j — 1 : length(set) 

* Signi = Deter _sign(f ,real_roots{t},nj,l); 

* If sign± = —1, then Index = Index U {s}, break; 

* End If 

— End for 

• End if 

3. End for 
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4. real^roots = real-roots \ real -roots(Index) . 



The following algorithm is to remove the intervals which satisfy p in (4). 
Algorithm 5: Dele-Posi 

Input: Two polynomial systems: / and p; isolated intervals real -roots 
of/; 

Output: real-roots = {X G real-roots\p(x) > 0,f(x) = 0,where& G 

X}. 

1. real-roots = Dele-Inequ(f,real_roots, p); 

2. real-roots = Dele-N onnega{f, real-roots, p); 

Until now, we have described all the parts of the algorithm for isolating 
the real roots of semi-algebraic system. Here, we present the whole algorithm 
in the following. 

Algorithm 6: real -root- semi 

Input: A semi- algebraic system in the form of (4). 
Output: isolated real root intervals of (4) or return {}; 

1. real-roots = real-root-isolate(f); 

2. If length(real-roots) = 0, then return {}; 

3. Else 

• real-roots = Dele-Inequ(f , real-roots, h); 

• If length(real -roots) = 0, then return {}; 

• Else 

— real-roots = Dele-N onnega(f , real-roots, n); 

— If length(real -roots) = 0, then return {}; 

- Else 

* real-roots = Dele-Posi(f , real-roots, p); 

* If length{r-eal -roots) = 0, then return {}; 

* End If 

- End If 

• End If 

4. End If 

5. Return real-roots. 
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In the following, a simple example is provided to illustrate algorithm 6 
for isolating the real roots of a semi-algebraic system. 

Example 3.2. Given a semi- algebraic system, 



Step 1: Denote the isolated real root intervals of [fi, f 2 ] by X, using 
algorithm 1, we obtain six intervals: 

X l = [3.419999999999901,3.420000000000101], 
[3.202328957744943,3.202328957745143]; 

X 2 = [3.419999999999901,3.420000000000101], 

[ - 17.753578957745148, -17.753578957744949]; 
X 3 = [ - 0.099812907756173, -0.099812907755973], 

[0.857641152899082,0.857641152899282]; 
X 4 = [ - 8.128471753948865, -8.128471753948666], 

[ - 7.758843802156950, -7.758843802156750]; 
X 5 = [0.584413218566582,0.584413218566782], 

[ - 2.373999555338525, -2.373999555338325]; 
X 6 = [0.560588744512268,0.560588744512469], 

[0.376338245290841, 0.376338245291041] . 

Step 2: Substituting X into x 1 — 3.42, we have 
eVl = [-0.000000000000099,0.000000000000100], 
ev 2 = [-0.000000000000099,0.000000000000101], 
ev 3 = [-3.519812907756173,-3.519812907755973], 
ev A = [-11.548471753948863,-11.548471753948665], 
ev 5 = [-2.835586781433418,-2.835586781433218], 
ev 6 = [-2.859411255487732,-2.859411255487531]. 

It is easy to see that x\ — 3.42 ^ at any point contained in intervals 
X 3 , X4, X§, Xq, so these intervals should be retained. Using algorithm 3 to 
process intervals X 1 and X 2 , we know that x\ — 3.42 is equal to zero at the 
solutions of {/1, / 2 }, so we discard Xi and X 2 at this step. 




h 
h 



= 5 + 13xi - 10x 2 - 82xj + 71xix 2 + 16x1 
= 403.22xx - 314.64 + 73.16a;? - 269.26x 1 a; 
48x? + 53x^ 2 - 2Sx^x 1 + 95.76a;2. 
x 1 - 0.3 > 0,x 2 > 0, Xl - 3.42 ^ 0. 
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Step 3: Substituting X 3 X 4 X 5 and X 6 into P = {xi — 0.3, x 2 }, we obtain 
the following four intervals 

P(X 3 ) = [ - 0.399812907756173, -0.399812907755973], 
[0.857641152899082,0.857641152899282]; 

P(X 4 ) = [ - 8.428471753948864, -8.428471753948665], 
[ - 7.758843802156949, -7.758843802156748]; 

P(X 5 ) = [0.284413218566582,0.284413218566782], 

[ - 2.373999555338525, -2.373999555338325]; 

P{X 6 ) = [0.260588744512268,0.260588744512469], 
[0.376338245290841,0.376338245291041]. 

where P(Xj) means interval extension P at X^, i = 3, 4, 5, 6. Obviously, X 3 , 
X 4 and X 5 should be deleted from X. 

Therefore, the isolated real root interval of SAS is 

X 6 = [0.560588744512268,0.560588744512469], 
[0.376338245290841, 0.376338245291041] . 

3.3. Comparison experiments 

Algorithm 6 has been implemented in Matlab 2010a named real_roots_semi. 
In this subsection, we do some experiments (see Appendix A) to compare it 
with Maple package Discover described in [TB] on the platform of Maple 15 
classic work sheet. All experiments are carried out on a Dell PC with Intel 
Core i3-2120 at 3.30GHZ and 4GB of RAM in Windows 7. 

Table 1 gives the comparison of execution time. In the first row, M_d, 
N_v and N_p denote the max degree, number of variables, number of real 
points of the semi-algebraic system. And oo means the program out of mem- 
ory. Among these examples, A. 5 is from [17], the other examples are chosen 
randomly. As shown in Table 1, Discover is faster than real^root^semi for 
Example A.l and A. 5, since these semi- algebraic system are simple. When 
the system becomes more complicated, real_root_semi is more efficient than 
Discover. From table 1 we can also find that the costs of example A. 6, A. 7 
and A. 9 are much higher than those of other examples compared because 
these examples can satisfy the inequations. 
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4. Applications 

Most engineering problems can be reduced to solving nonlinear equations. 
In this section, some applications from robot kinematic, piecewise algebraic 
curve, chemical engineering, circuit design are tested to show the efficiency 
of our algorithm. 

4-1. Polynomial system 

In this subsection, we first investigate applications of our algorithm to 
some polynomial systems. 

Example 4.1. Production of synthesis gas in an adiabatic reacto r[34\[ . 

Xix-j + 2x 2 x 7 + x 3 x 7 — 2x 6 = 0. 
x 7 (x 3 + x 4 + 2x + 5) - 2 = 0. 
7xi + 7x2 + 7^5 — 1 = 0. 
xi + x 2 + x 3 + x 4 + x 5 - 1 = 0. 

400x^1 - 178370x3X5 = 0. , . 

' X1X3 - 2.60582 2 X4 = 0. ^ > 

-28837xix 7 - 139009x 2 x 7 - 78213x 3 x 7 + 18927x 4 x 7 + 8427x 5 x 7 + 
13492 - 10690x 6 = 0. 
< Xi < l,i = 1,2,3,4,5. 
< Xi < 5, % = 6, 7. 

The above system of equations represents three atom balances, a mole 
fraction constraint, two equilibrium relations, and an energy balance equa- 
tions. We obtain 8 real roots of this equation through 0.536665 seconds 



Example No. 


Max degree 


N_v 




Discover 


real-roots _ semi 


A.l 


2 


3 


1 


0.062 


0.259715 


A.2 


3 


3 


1 


0.483 


0.188181 


A.3 


4 


3 





219.759 


0.461846 


A.4 


2 


4 


1 


13.947 


0.302547 


A.5 


3 


4 


1 


0.31 


0.756411 


A.6 


4 


4 


3 


9.734 


6.106957 


A.7 


2 


4 


2 


13.821 


1.942310 


A.8 


4 


5 


5 


00 


0.794765 


A.9 


2 


6 


5 


00 


11.272852 



Table 1: Time comparison, units:s. 
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and 8 interval iterations for the 8 solution, then remove 7 real roots after 
0.061837 seconds. Let the midpoint of interval be the approximate root of 
the problem, then 

x = [0.322870839476541,0.009223543539188,0.046017090960632, 
0.618171675070824, 0.003716850952815, 0.576715395935549, 2.977863450791145]. 

Example 4.2. Robot kinematic problem [35]. 

' (0.004731x3 - 0.1238)xi - (0.3578x 3 + 0.001637)x 2 + x 7 - 0.9338x 4 
-0.3571 = 0. 

0.2238xix 3 + 0.7623x 2 x 3 + 0.2638X! - x 7 - 0.07745x 2 - 0.6734x 4 
-0.6022 = 0. 

x 6 x 8 + 0.3578 + 0.004731x 2 = 0. 
< -0.7623xi + 0.2238x2 + 0.3461 = 0. (7) 
x\ + x\ - 1 = 0. 
x\ + x\ - 1 = 0. 
x\ + x\ - 1 = 0. 

x 7 + x i ~ i- 

k Xi > -l,Xi < l,i = 1,2, • • • ,8. 

All 16 real solutions were founded in 0.826805 seconds of CPU time, which 
is more efficient than 3.783 seconds CPU time in [37] to solve the problem. 

Example 4.3. This example addresses the equilibrium of the products of 
a hydrocarbon combustion process [36]. The problem is reformulated in the 
1 element variables' space. 

ym + vi - 3^/5 = o. 

2yi2/ 2 + yi + 3i2ioi/l + 2/22/I + R7V2V3 + R9V2V4 + RsV2 - Ry5 

22/22/1 + #72/22/3 + 2i?52/3 2 + #62/3 - % 5 = 0. 

< #92/22/4 + 2yj - ARy 5 = 0. 

2/12/2 + 2/1 + #102/1 + 2/22/1 + #72/22/3 + #92/22/4 + #82/2 + #52/3 

+#62/3 + 2/1-1=0. 
k Vi - 0.0001 > 0, 100 - yi > 0, % = 1, 2, • • • , 5. 

The value of parameter Ri,i = 1,2, ■■ ■ ,10 can be found in [36]. The 
method in [37J finds the single solution after 217.7 seconds of CPU time. 
Using our method, we only use 0.2964 seconds CPU time. 

Example 4.4. This example for computing the real roots of piecewise alge- 
braic curve is modified from /2J/. 
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-1 



Fig. 1 . Two piecewise algebraic curves 



Let A = {a, b, c, d} be a regular triangulation of rectangular domain 
ABCD in R 2 , where a = [AOB], b = [BOC], c = [COD] and d = [AOD], 
where A = (-3, 0), B = (0, -3), C = (3, 0) D = (0, 3) and O = (0,0), f and 
g are algebraic curves (See Fig. 1). 

Suppose that f,g £ S\ (A) and 



We take the cell a for instance, isolating the intersection points of {f\\ a , gi\ a } 
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lying in a can be converted to isolating the real roots of the following semi- 
algebraic system: 

( x 3 — y 3 + 3y 2 + xy — 3x — 4 = 0. 



[ -x > 0, -y > 0,x + y + 3 > 

Using the algorithm real -root- semi, we know that there is only one real root 
lying in a. Similarly, we conclude that {/i| c , <?i| c } has one real point in the 
interior of c, and respectively, {/i|&, <7i|&} and {fi\d, 9i\d} has two common 
points in the cell of b and d. 

Let the polynomial {fi\b, <?i|&} be inequations, and append it to system 
(9), we conclude that a and b have a common real point. Hence, / and 
g have 5 real points. Methods given by Wang, et al., based on interval 
method, Groebner basis in [H [21 El IH E] can also find the commom points of 
piecewise algebraic curve, but they did not consider these points which lie in 
the boundary. 

4-2. Transcendental functions system 

In practice, systems of transcendental functions appear in many applica- 
tions [3H E5], e.g., logarithm functions, exponential functions and trigono- 
metric functions. 

As far as we know, there exists three kinds of numerical methods to solve 
nonlinear equations: homotopy continuation method, Newton's iteration and 
interval bisection method. Each kind of algorithm has its advantages and 
disadvantages. Homotopy continuation method is a global convergence al- 
gorithm, but it can only be used for solving polynomial equations. Given 
an appropriate initial point, the convergence of Newton's method is typical 
local quadratic, but how to choose the initial point is a difficult problem. 
The interval method is unacceptable in computation when some conditions 
are complicated, such as more variables, the huge width of the initial interval 
and the more accuracy of the final result. 

The basic idea of our method is to replace the exponential function by 
its Taylor expansion first. In this step, rather than using the initial interval, 
we reduce the width of the initial interval using interval algorithm, until the 
errors between the exponential function and its Taylor expansion is smaller 
than a tolerance. We then apply homotopy method to the polynomial system 
to obtain the approximate real roots. At last, the real roots of transcendental 




(9) 
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functions are found based on these approximate real roots using the Newton's 
method. 

In the following, we apply this idea to solve a circuit design problem from 
[38]. 

Example 4.5. Circuit design problem with extraordinary sensitivities to 
small perturbations. 

' (1 - XlX2 ) X3 (e(^^"-^^ 10_3 -f 5 ^8io- 3 )) - l) - g 5k + g AkX2 = o. 

(1 - x^x^e^ 9 ^- 9 ^- 9 ^ 10 '^ 9 ^ 910 ' 3 ^ - 1) - g 5k xi + 9ik = 0. 
Xix 3 — x 2 x$ = 0. (10) 
k = 1,2,3,4. 
. Xi E [0,10],z = I,-- - ,9. 

The constant g^. are given by the following matrix 

/ 0.4850 0.7520 0.8690 0.9820 \ 

0.3690 1.2540 0.7030 1.4550 

5.2095 10.0677 22.9274 20.2153 

23.3037 101.7790 111.4610 191.2670 

\ 28.5132 111.8467 134.3884 211.4823 / 

Step 1: Since the exponential function in system (10) only contains vari- 
ables Xi,i = 5,6,7,8,9, we only bisect the intervals of these variables by 
interval Newton iteration. This is different from other algorithms. After 10 
times interval iterations for each variable, we have the last interval of these 
variables 

[7.998046875,8.007812500], 
[7.998046875,8.007812500], 
[5.000000000,5.009765625], 
[0.996093750,1.005859375], 
[1.992187500,2.001953125]. 

So we only need 50 bisections, and this progress takes about 9 seconds. 
Now, we compute the interval extension of git = x$(gik — gsfcXylO -3 — 
5f 5fc x 8 10~ 3 ) and g 2k = x 6 (g lk - g 2k ~ g3 k x 7 W~ 3 + fi^glO" 3 ), and using the 
Taylor expansion of the exponential function in the midpoint of the interval 
extension of g\ k and g 2k to replace e 9lk and e 92k for k = 1, 2, 3,4 in equation 
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(10). At this point, we use the second order Taylor expansion for e 914 and the 
first order for other exponential functions. Note that the errors generated by 
this step are not more than 1CT 2 . 

Step 2: For the polynomial system obtained from step 1, there exists 37 
homotopy paths, but only 4 roots are founded using homotopy algorithm, 
and only 2 of them are real roots. This process takes 1.469615 seconds. 

Step 3: Let the two real points be the initial point of the Newton's it- 
eration for system (10). After 0.092 seconds we find that this two points 
converge to the same point. 

x = [.899999952618, .449987471886, 1.00000648241, 7.99997144063, 5.00003127610, 
.999987723423, 2.00006854190, 7.99969268290, 2.00005248366]. 

Compared with x, we find that the errors between x and the initial point 
from step 2 is 1.2602 and 0.0017, respectively. The last one is close to x 
which ensure the convergence of the Newton iterative algorithm in this step. 

From the above description, we totally use about 11 seconds to find the 
real point of equation (10), which is much less than 436.5 seconds in [37]. 

5. Conclusion 

For a class of semi-algebraic system, this paper presents a numerical 
method for isolating the real roots. The algorithm first obtains the isolated 
intervals of a zero-dimensional polynomial system using hybrid technique 
[2"2"] . and substitute it into the inequations, nonnegative system, and positive 
systems to remove some intervals. Such implementation will be efficient if 
the interval extension of constrained equation does not contain zero at the 
isolated real root interval of the system. Otherwise, we give a complete nu- 
merical algorithm to determine if a polynomial equation is equal to zero. 
We implement our algorithm in Matlab environment. Many random exam- 
ples have been checked along with comparison to famous software Discover 
. At last, this new method is used to solve problems originated from several 
applications. For transcendental functions, we give an idea for solving this 
problem which is a combination of Newton's method, homotopy algorithm 
and interval arithmetic. However, the number of bisections, the order of 
Taylor expansion and the error control are deserved to further study. 
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Appendix A. 



-3 + 95x 3 + 29xix 2 + 63x1X3 + 37x^ - 3Qx 2 3 = 
— 2xi + 46x 2 + 4xix 2 + 50xix 3 + 61x 3 x 2 — 99xg : 
72 - 49xx + 13x 2 - 15x? - 12x5; + 23X 2 = 0. 
3x 2 + 2xi - 1 > 0, 3xi - 4x 2 - 3 > 0. 

' -85xi - 76x 2 + 75x 3 - 41x 2 - 84x| = 0. 

-39 + 26xix 3 + 47x? - 47x 2 x 2 + 91x^xi + 43x 3 xix 2 = 0. 

16xi + 25x 3 + 13xix 2 + x\ — 8x 3 x 2 — 74x 2 x^ = 
< 19x 3 + 10xix 3 + 53x 3 x 2 - 97x 3 x 2 + 57x^x x + Q8x 2 x 2 3 > (A.2) 

47xiX2 + 58x3X 2 + 30x| + 55x3xf + 56x3X1X2 + 55x| > 

48x 2 + 34x 2 + 32xf - 41x|xi - 5Qx 3 2 + Q3x 2 x 2 3 > 
k 27 + llxi + 46x 3 - 42xix 3 - 60x? - 48x|x 3 > 

-49x 2 + 31x 3 x 2 + 73x1X3 + 95x3X1 + 68xix 2 x^ - 29x1x3 s = 0. 
37 + 5x? _ 3Qx 2 _ 57x2X j + g 5xiX s + 80x 2 x 2 = 0. 

30x 2 — 3x 3 x 2 — 56xfx 2 — 91xlxl — 70x 3 x 2 xi + 42x| = 
< 62 + 96x 2 - 51x 3 + 89x 2 + 14x^xi - 79xlx 3 > (A.3) 
86 + 57xi - 35x 2 + 57X 2 , + 28x^ + 63x 3 x 2 > 
-61xi + 45x 3 - 50x 2 2 - 22x1 + 48x 2 x 2 - 82x 2 xj > 
-85x 3 x 2 - 44x 3 x 2 - 31x^X1 + 45xix 3 x 2 + 49x| - 58x^ > 



(Al) 
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94 - 99xix 3 + 46x 2 x 4 + 41xj + 95x 3 x 4 = 0. 

— 59xi — 15x4 + 30xix 2 — 57xix 3 + 50x1 + 62x2^3 = 0. 

53 - 9xi - 99x 4 + 12xix 4 + 20x 2 x 4 + 99xf = (A.4) 

-82 - 22x 3 - 60x^ + 3x 2 x 4 - 83x 3 x 4 - 84x| = 

-74x 2 x 3 + 87x 3 + 68x 4 x^ - 88x 3 x^ - 88x 2 xj - 87x 2 x 3 x 4 > 



2 — 6.5xi + x^x 2 — 0.5x 3 = 0. 
6x1 — x\x2 — 5x4 + 5x 2 = 0. 
< 2 - 6.5x 3 + x§x 4 - 0.5xi = (A.5) 
6x 3 — x 3 X4 + 1 + 0.5x 2 — 0.5x4 = 
Xi > 0, x 2 > 0, x 3 > 0, x 4 > 



9x^ — 20 + x^x 4 — 5x4 — x\ — Ax\ x 2 + 20xix 2 — 4x^x 2 x 4 + 20x 2 X4 
+3xjxl - 15x% = 0. 

3xfx 2 — 39x5 + 4x4X3 ~~ 52x x x 2 — 5x 3 XiX 2 + 65x]X 3 + 3xix|x 4 

— 7 — 7x 2 + 4xix 3 — 2x4X4 + 5x 2 x 3 + 5x 2 X4 = 0. 
3 + 2xi + 3x 3 — 3x^ + 5xix 2 — x 3 X4 = 0. 
x\ - 5 ^ 0, x\ - 13 ^ 0. 

(A.6) 



2x5 ~~ 3x^x 2 — 8x5 — 5x1X3 — 12xix 2 + 5xfx 3 — 12xix 3 x 2 — 20xix 3 
+4x^X4 — 12xiX 2 X4 — I6X1X4 — 12x| — I6X2 — 9x 2 x 3 — 12x 2 x 3 = 0. 
50x 3 — 65 — 20xix 3 + 26xi — 30xix 3 x 2 + 39xix 2 — 70x 2 x| + 91x 2 x 3 
-70a:| + 91x| - 40x|x 4 + 52x 3 x 4 = 0. 
5xi — 2x 3 — 2xiX 4 — 2x 2 x 4 — x 3 — 7x^ = 0. 
—7 — 7xi + 2xix 2 — 7x 2 x 4 — 3x 3 x 4 + 5x| = 0. 
xi - 3x 2 -4 ^ 0,x 3 - 1.3 ^ 0. 

(A.7) 
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66x4X 2 xi — 44x1x2X4 + 66x1X20:4 — I65X4X5X2 + 110x4x|xi — I65X4X5 = 0. 

-10xi + 62x 3 - 82x 4 + 80x 5 - Ux 2 2 + 71x 3 x 4 = 0. 

74xiX 2 + 72xiX 4 + 37x 2 x 4 — 23x 2 x 5 + 87x 3 x 5 + 44x^ = 0- 

11 - 49x 2 - 47x 5 + 40xj - 81x 2 x 3 + 91x 2 x 4 = 0. 

— 28xi + 16x 2 + 30x4 — 27xix 2 — 15xix 4 — 59x 3 X4 = 0. 

3x2-2xi + 3^0,x 4 ^0. 

(A.8) 



36x 2 + 91x 3 — 22x 2 x 3 + 51x 3 x 6 — 27x 4 x 5 + 50x1? = 0. 
25xi + 31x 2 — 27x 4 + 65x]X 2 + 88x 2 x 5 + 10x 3 x 5 = 0. 
95x 2 + 68xiX 2 — 29xxx 6 + 5x 2 x 4 — 26xg — 51x 3 x 4 = 0. 
5 — 36x4 — 57xix 2 + 85x 3 + 8OX4X6 + 90xg = 0. 
65xix 2 — 12xix 5 + 78xiX6 + 5x 2 — 63x 2 X6 — 5x 3 xs = 0. 
-70 + 42xx + 9x1 - 21xxx 2 - 27x x x 5 - 79x^ = 0. 
x\ + x 5 x 6 — 2 > 0. 
3xi - 4x 2 + 4 > 0. 
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